{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#My name: Joe Puccio\n", "#Collaborators: Fred Landis, Alan Wu, Zen Yang\n", "\n", "#Helper functions courtesy of Master McMillan\n", "import string\n", "import array as arrayModule\n", "import sys\n", "from numpy import *\n", "import time \n", "\n", "# The following classes are used to construct \n", "# and represent Suffix trees\n", "\n", "class Edge:\n", " count = 0\n", " def __init__(self, dstNode, first, last):\n", " self.dstNode = dstNode\n", " self.first = first\n", " self.last = last\n", " Edge.count += 1\n", " def split(self, suffix, suffixTree):\n", " # add new explicit node\n", " newIndex = len(suffixTree.nodes)\n", " suffixTree.nodes.append(0)\n", " # add new suffix edge\n", " newFirst = self.first + len(suffix)\n", " suffixTree.edgeLookup[newIndex, suffixTree.string[newFirst]] = Edge(self.dstNode, newFirst, self.last)\n", " # shorten this edge\n", " self.last = newFirst - 1\n", " self.dstNode = newIndex\n", " return newIndex\n", " def isLeafEdge(self, suffixTree):\n", " return (self.last == suffixTree.lastIndex)\n", " def __len__(self):\n", " return self.last - self.first + 1\n", " def __repr__(self):\n", " return \"Edge(%d, %d, %d)\" % (self.dstNode, self.first, self.last)\n", "\n", "class Suffix:\n", " def __init__(self, srcNode, first, last):\n", " self.srcNode = srcNode\n", " self.first = first\n", " self.last = last\n", " def is_explicit(self):\n", " return self.first > self.last\n", " def is_implicit(self):\n", " return self.first <= self.last\n", " def canonicalize(self, suffixTree):\n", " if self.is_implicit():\n", " edge = suffixTree.edgeLookup[self.srcNode, suffixTree.string[self.first]]\n", " if (len(edge) <= len(self)):\n", " self.first += len(edge)\n", " self.srcNode = edge.dstNode\n", " self.canonicalize(suffixTree)\n", " def __len__(self):\n", " return self.last - self.first + 1\n", " def __repr__(self):\n", " return \"Suffix(%d, %d, %d)\" % (self.srcNode, self.first, self.last)\n", "\n", "class SuffixTree:\n", " def __init__(self, string):\n", " self.string = string # save a pointer to the string\n", " self.alphabet = set() # alphabet of string\n", " self.nodes = arrayModule.array('l') # index of ith node's parent\n", " self.nodes.append(0) # add root node\n", " self.edgeLookup = {} # adjacency list indexed by (srcNode, char)\n", " self.lastIndex = len(string) - 1\n", " activePoint = Suffix(0, 0, -1)\n", " for i in xrange(len(string)):\n", " self.alphabet.add(string[i])\n", " self.addPrefix(i, activePoint)\n", " def addPrefix(self, last, activePoint):\n", " LastParentNode = -1\n", " while True:\n", " ParentNode = activePoint.srcNode\n", " if activePoint.is_explicit():\n", " if (activePoint.srcNode, self.string[last]) in self.edgeLookup:\n", " break\n", " else: #potentially split an implicit node\n", " edge = self.edgeLookup[activePoint.srcNode, self.string[activePoint.first]]\n", " if (self.string[edge.first + len(activePoint)] == self.string[last]):\n", " break\n", " else:\n", " ParentNode = edge.split(activePoint, self)\n", " self.nodes.append(-1)\n", " self.edgeLookup[ParentNode, self.string[last]] = Edge(len(self.nodes)-1, last, self.lastIndex)\n", " # add suffix link\n", " if (LastParentNode > 0):\n", " self.nodes[LastParentNode] = ParentNode\n", " LastParentNode = ParentNode\n", " if (activePoint.srcNode == 0):\n", " activePoint.first += 1\n", " else:\n", " activePoint.srcNode = self.nodes[activePoint.srcNode]\n", " activePoint.canonicalize(self)\n", " if (LastParentNode > 0):\n", " self.nodes[LastParentNode] = ParentNode\n", " activePoint.last += 1\n", " activePoint.canonicalize(self)\n", " def leafIndices(self, nodeIndex=0, lenSoFar=0):\n", " indexList = []\n", " if (self.nodes[nodeIndex] < 0):\n", " indexList.append(self.lastIndex + 1 - lenSoFar)\n", " else:\n", " for char in self.alphabet:\n", " try:\n", " edge = self.edgeLookup[nodeIndex, char]\n", " if edge.isLeafEdge(self):\n", " indexList.append(self.lastIndex + 1 - len(edge) - lenSoFar)\n", " else:\n", " indexList += self.leafIndices(edge.dstNode, lenSoFar + len(edge))\n", " except KeyError:\n", " continue\n", " return indexList\n", " def distinct(self, nodeIndex=0, lenSoFar=0):\n", " distinctList = []\n", " # examine children of node\n", " for char in self.alphabet:\n", " try:\n", " edge = self.edgeLookup[nodeIndex, char]\n", " if edge.isLeafEdge(self):\n", " distinctList.append((edge.first-lenSoFar, lenSoFar+1))\n", " else:\n", " distinctList += self.distinct(edge.dstNode, lenSoFar + len(edge))\n", " except:\n", " continue\n", " return distinctList\n", " def thread(self, target):\n", " nodeIndex = 0 # start from root\n", " i = 0 # characters threaded so far\n", " while (i < len(target)):\n", " try:\n", " edge = self.edgeLookup[nodeIndex, target[i]]\n", " prefix = self.string[edge.first:edge.last+1]\n", " if (target[i:].startswith(prefix) or prefix.startswith(target[i:])):\n", " i += len(prefix)\n", " else:\n", " return []\n", " nodeIndex = edge.dstNode\n", " except KeyError:\n", " return []\n", " return self.leafIndices(nodeIndex, i)\n", " def printTree(self, nodeIndex=0, prefixLength=0, prefix=\"\"):\n", " branches = list()\n", " for c in self.alphabet:\n", " try:\n", " edge = self.edgeLookup[nodeIndex, c]\n", " extent = prefix + \"(%d)-%s-\" % (nodeIndex, self.string[edge.first:edge.last+1])\n", " subtree = self.printTree(edge.dstNode, prefixLength+len(edge), extent)\n", " if (len(subtree) > 0):\n", " branches.append(subtree)\n", " elif (edge.isLeafEdge(self)):\n", " print extent + \"(%d)[%d]\" % (edge.dstNode, edge.first-prefixLength)\n", " except KeyError:\n", " pass\n", " return branches\n", " \n", "# a one-line function to compute a Suffix array\n", "def argsort(text):\n", " return sorted(range(len(text)), cmp=lambda i,j: -1 if text[i:] < text[j:] else 1)\n", "\n", "def findFirst(pattern, text, sfa):\n", " \"\"\" Finds the index of the first occurence of pattern in the suffix array \"\"\"\n", " hi = len(text)\n", " lo = 0\n", " while (lo < hi):\n", " mid = (lo+hi)//2\n", " if (pattern > text[sfa[mid]:]):\n", " lo = mid + 1\n", " else:\n", " hi = mid\n", " return lo\n", "\n", "def findLast(pattern, text, sfa):\n", " \"\"\" Finds the index of the last occurence of pattern in the suffix array \"\"\"\n", " hi = len(text)\n", " lo = 0\n", " m = len(pattern)\n", " while (lo < hi):\n", " mid = (lo+hi)//2\n", " i = sfa[mid]\n", " if (pattern >= text[i:i+m]):\n", " lo = mid + 1\n", " else:\n", " hi = mid\n", " return lo" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "chr1 = open(\"HumChr01.seq\",'rb').read()\n", "import sys\n", "sys.setrecursionlimit(10000)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 100, 10, 0.001 secs 0, 0.000 secs\n", " 1000, 1855, 0.026 secs 3, 0.000 secs\n", " 10000, 16607, 0.218 secs 112, 0.001 secs\n", " 100000, 165240, 1.916 secs 1823, 0.011 secs\n", " 1000000, Crashed \n", "\n", "Number of CATs found: 1823\n" ] } ], "source": [ "#Problem 1\n", "\n", "N = 100\n", "treeTime = 0.0\n", "while (N <= 100000000) and (treeTime < 100.0):\n", " print \"%10d, \" % N,\n", " text = chr1[10000:10000+N]\n", " start = time.clock()\n", " try:\n", " sTree= SuffixTree(text)\n", " except RuntimeError:\n", " print \"Crashed \\n\"\n", " print \"Number of CATs found:\", len(catsfound)\n", " break\n", " treeTime = time.clock() - start\n", " print \"%8d, %6.3f secs\" % (len(sTree.nodes), treeTime),\n", " start = time.clock()\n", " catsfound = sTree.thread(\"CAT\")\n", " findTime = time.clock() - start\n", " print \"%8d, %6.3f secs\" % (len(catsfound), findTime)\n", " N*=10\n", " " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 100, 100, 0.000 secs 0, 0.000 secs\n", " 1000, 1000, 0.009 secs 3, 0.000 secs\n", " 10000, 10000, 0.167 secs 112, 0.000 secs\n", " 100000, 100000, 9.037 secs 1823, 0.000 secs\n", " 1000000, 1000000, 997.981 secs 13849, 0.001 secs\n", "\n", "Number of CATs found: 13849\n" ] } ], "source": [ "#Problem 2\n", "N = 100\n", "arrayTime = 0.0\n", "while (N <= 100000000) and (arrayTime < 100.0):\n", " print \"%10d, \" % N,\n", " text = chr1[10000:10000+N]\n", " start = time.clock()\n", " sArray = argsort(text)\n", " arrayTime = time.clock() - start\n", " print \"%8d, %6.3f secs\" % (len(sArray), arrayTime),\n", " start = time.clock()\n", " lo = findFirst(\"CAT\", text, sArray)\n", " hi = findLast(\"CAT\", text, sArray)\n", " findTime = time.clock() - start\n", " print \"%8d, %6.3f secs\" % (hi-lo, findTime)\n", " N *=10\n", "print \"\\nNumber of CATs found:\", hi-lo" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 10000, 0.134 secs 0.342 secs\n", " 11000, 0.129 secs 0.422 secs\n", " 12000, 0.121 secs 0.246 secs\n", " 13000, 0.141 secs 0.354 secs\n", " 14000, 0.165 secs 0.392 secs\n", " 15000, 0.234 secs 0.342 secs\n", " 16000, 0.198 secs 0.588 secs\n", " 17000, 0.218 secs 0.427 secs\n", " 18000, 0.256 secs 0.439 secs\n", " 19000, 0.283 secs 0.541 secs\n", " 20000, 0.312 secs 0.503 secs\n", " 21000, 0.396 secs 0.491 secs\n", " 22000, 0.385 secs 0.543 secs\n", " 23000, 0.396 secs 0.540 secs\n", " 24000, 0.451 secs 0.714 secs\n", " 25000, 0.465 secs 0.595 secs\n", " 26000, 0.509 secs 0.606 secs\n", " 27000, 0.568 secs 0.706 secs\n", " 28000, 0.609 secs 0.638 secs\n", " 29000, 0.704 secs 0.753 secs\n", " 30000, 0.742 secs 0.666 secs\n", " 31000, 0.769 secs 0.786 secs\n", " 32000, 0.755 secs 0.846 secs\n", " 33000, 0.802 secs 0.759 secs\n", " 34000, 0.919 secs 0.863 secs\n", " 35000, 0.985 secs 0.874 secs\n", " 36000, 1.066 secs 0.792 secs\n", " 37000, 1.104 secs 0.915 secs\n", " 38000, 1.035 secs 0.937 secs\n", " 39000, 1.088 secs 0.988 secs\n", " 40000, 1.307 secs 1.275 secs\n", " 41000, 1.332 secs 1.126 secs\n", " 42000, 1.393 secs 1.312 secs\n", " 43000, 1.460 secs 1.338 secs\n", " 44000, 1.588 secs 1.426 secs\n", " 45000, 1.625 secs 1.769 secs\n", " 46000, 1.740 secs 1.600 secs\n", " 47000, 1.797 secs 1.790 secs\n", " 48000, 1.870 secs 1.607 secs\n", " 49000, 1.809 secs 1.169 secs\n", " 50000, 1.836 secs 1.095 secs\n", " 51000, 1.871 secs 1.304 secs\n", " 52000, 2.075 secs 1.201 secs\n", " 53000, 2.098 secs 1.153 secs\n", " 54000, 2.246 secs 1.391 secs\n", " 55000, 2.308 secs 1.273 secs\n", " 56000, 2.354 secs 1.305 secs\n", " 57000, 2.339 secs 1.393 secs\n", " 58000, 2.635 secs " ] } ], "source": [ "#Problem 3\n", "import matplotlib.pyplot as plt\n", "\n", "class SuperTree(SuffixTree):\n", " def getSuffixArray(self, nodeIndex=0, prefixLength=0):\n", " if (nodeIndex == 0):\n", " self.suffixArray = []\n", " for c in sorted(self.alphabet):\n", " try:\n", " edge = self.edgeLookup[nodeIndex, c]\n", " if (edge.isLeafEdge(self)):\n", " self.suffixArray.append(edge.first-prefixLength)\n", " else:\n", " self.getSuffixArray(edge.dstNode, prefixLength+len(edge))\n", " except KeyError:\n", " pass\n", " return self.suffixArray\n", "\n", " \n", "N = 10000\n", "arrayTime = 0.0\n", "arrayTimeArray = []\n", "treeTimeArray = []\n", "NArray = []\n", "\n", "while(arrayTime < 2.4):\n", " NArray.append(N)\n", " print \"%10d, \" % N, \n", " text = chr1[10000:10000+N]\n", " start = time.clock()\n", " sArrray = argsort(text)\n", " arrayTime = time.clock() - start\n", " arrayTimeArray.append(arrayTime)\n", " print \"%6.3f secs\" % (arrayTime), \n", " start = time.clock()\n", " sTree = SuperTree(text)\n", " tArray = sTree.getSuffixArray()\n", " treeTime = time.clock() - start\n", " treeTimeArray.append(treeTime)\n", " print \"%6.3f secs\" % (treeTime)\n", " N+=1000\n", "\n", "plt.plot(NArray, arrayTimeArray, 'bs', NArray, treeTimeArray, 'g^')\n", "plt.show()\n", "#using this plot, we see that the length of a text where first constructing \n", "#a suffix tree will be faster than the given algorithm of suffix array construction\n", "#is around 32000 characters. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.5" } }, "nbformat": 4, "nbformat_minor": 0 }